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We propose a new numerical method to compute quasi-equilibrium sequences of general relativistic 
irrotational binary neutron star systems. It is a good approximation to assume that (1) the binary 
star system is irrotational, i.e. the vorticity of the flow field inside component stars vanishes every- 
where (irrotational flow), and (2) the binary star system is in quasi-equilibrium, for an inspiraling 
binary neutron star system just before the coalescence as a result of gravitational wave emission. 
We can introduce the velocity potential for such an irrotational flow field, which satisfies an elliptic 
partial differential equation (PDE) with a Neumann type boundary condition at the stellar surface. 
For a treatment of general relativistic gravity, we use the Wilson-Mathews formulation, which as- 
sumes conformal flatness for spatial components of metric. In this formulation, the basic equations 
are expressed by a system of elliptic PDEs. We have developed a method to solve these PDEs 
with appropriate boundary conditions. The method is based on the established prescription for 
computing equilibrium states of rapidly rotating axisymmetric neutron stars or Newtonian binary 
systems. We have checked the reliability of our new code by comparing our results with those of 
other computations available. We have also performed several convergence tests. By using this code, 
we have obtained quasi-equilibrium sequences of irrotational binary star systems with strong gravity 
as models for flnal states of real evolution of binary neutron star systems just before coalescence. 
Analysis of our quasi-equilibrium sequences of binary star systems shows that the systems may not 
suffer from dynamical instability of the orbital motion and that the maximum density does not 
increase as the binary separation decreases. 



I. INTRODUCTION 

At the beginning of the next century, advanced grav- 
itational wave detectors (LIGO/VIRGO/TAMA/GEO; 
101) will provide us with signals from coalescing binary 
neutron star systems. According to the type of com- 
pact objects involved (for example, black holes/neutron 
stars, isolated/binary objects, solar mass/super massive 
objects), various theoretical methods have been devel- 
oped or are under development to predict and interpret 
gravitational wave signals which will be observed. Com- 
pact binary systems of solar mass sizes are also promising 
sources of 7— ray bursts or r-process nuclei. These are 
important issues of recent astrophysics (see for example 
1^ and references therein). 

When we consider the final evolutionary stage of the 
binary neutron star (BNS) as a result of gravitational 
wave emission, it is necessary to solve the internal struc- 
ture of the neutron star (NS) at the late inspiraling stage 
as well as the early stage of merging, since the size of the 
star is comparable to the separation of the two neutron 
stars. In this case, (1) tidal deformation of a component 
star, (2) the general relativistic effect on the orbital mo- 
tion, and (3) the general relativistic effect on the internal 



structure of a star are coupled to affect the physics of the 
binary system. Such finite size effects are not negligible 
for the last A/" ^ 10 cycles of the inspiraling phase just 
before coalescence 

According to the scenario of the final evolution of the 
close BNS ||], the flow field inside component stars just 
prior to coalescence becomes irrotational, i.e. a vorticity 
seen from the non-rotating observer's frame vanishes just 
prior to coalescence. This irrotational flow is realized 
because (1) the viscosity of neutron star matter is too 
weak to synchronize the spins of the component stars and 
the orbital angular velocity of the binary system during 
the evolution, and (2) the initial spin period of each NS 
is much larger than the final orbital period of the BNS. 
It is on precisely this stage of the binary evolution that 
we focus in the present paper. 

Full 3D numerical simulations are most promising to 
investigate such final phases of BNS systems. Numer- 
ical methods for fully dynamical 3D GR computations 
of BNS's are developing rapidly (see for example [§-^). 
However it is still difficult to perform simulations over 
the whole A/" ^ 10 cycles of the inspiraling BNS because 
there remain some technical problems to overcome for 
newly developed 3D GR codes, such as, maintaining sta- 
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bility and accuracy of numerical schemes and restrictions 
on the computational power of present computers. Re- 
cently, Shibata has succeeded in developing a new nu- 
merical method to compute dynamical evolution of the 
space-time and compact object (s) such as a single NS or 
a BNS system in a stable manner Q. Therefore simu- 
lations over several orbital periods of BNS will be per- 
formed in the near future. However, even after that, it 
is desirable to have totally different and supplementary 
methods to investigate such systems to test the reliability 
of the theory based on numerical methods. 

One of promising alternative ways to tackle this prob- 
lem is to compute quasi-equilibrium sequences of BNS 
configurations. Since the time scale of evolution as a re- 
sult of gravitational wave (GW) emission is longer than 
the time scale of the orbital motion even at this stage, 
quasi-equilibrium is a good approximation for a realistic 
BNS system Q]. Therefore, a quasi-equilibrium sequence 
with a constant rest mass constructed by changing the 
binary separation could be considered as an evolution- 
ary sequence of such an inspiraling BNS system. Also, 
quasi-equilibrium configurations will give accurate initial 
conditions for the future 3D GR dynamical computa- 
tions. Therefore, it is important to develop a numerical 
method to compute quasi-equilibrium configurations of 
close BNS's. 

In this paper we introduce a new numerical method to 
compute quasi-equilibrium configurations of irrotational 
BNS systems. It is a straightforward extension of the 
computational method used for the Newtonian irrota- 
tional binary systems developed by the present authors 
1^ . We discuss the formulation of the problem separately 
for (1) the gravitation part (the Einstein equations) and 
(2) the fluid part (the NS, i.e., the equation of motion, 
the continuity equation, the equation of state (EOS) and 
so on) . As the first step of the extension to general rela- 
tivistic strong gravity, a simplified system of the Einstein 
equations is employed. We use the formulation developed 
by Wilson and Mathews j|] and Baumgarte et al. jl^ , in 
which the metric is decomposed into (3-1-1) form and its 
spatial part is assumed to be conformally fiat. For the 
fiuid part, we solve the equations for irrotational flow 
whose formulation in GR has been developed by Shibata 
and Teukolsky Q (see also (l3)). O ur new com- 
putational code reproduces results of several indepen- 
dent computations previously made for GR co-rotating 
BNS's and Newtonian irrotational binary systems with 
reasonable accuracy. The convergence test of the numer- 
ical scheme is shown as well. Using this new code, we 
have constructed quasi-equilibrium sequences of irrota- 
tional BNS's with constant rest mass, which approximate 
the final evolutionary track of a realistic BNS as a result 
of GW emission. Analysis of our quasi-equilibrium se- 
quences of binary star systems suggests that the systems 
may not suffer from dynamical instability of the orbital 
motion and that the maximum density does not increase 
as the binary separation decreases. 

This paper is organized as follows. In section II we 



review the formulation for quasi-equilibrium irrotational 
binary systems in GR and derive the basic equations. 
In section III, we present the new numerical method for 
solving them. In section IV, we present several com- 
parisons with other computations, a convergence test of 
the new numerical scheme and the results for the binary 
configurations and quasi-equilibrium sequences. In sec- 
tion V, we discuss the validity of the present results and 
future prospects. 

II. REVIEW OF THE FORMULATION OF THE 
PROBLEM 

A. Quasi-equilibrium irrotational binary systems in 
GR: Summary of the assumptions and of the 
formulation 

Because of the finite size effect of stars in a close binary 
system, it is necessary to solve for the structures of the 
stars without approximation for the last J\f <, 10 cycles 
of the inspiraling phase just before coalescence. In gen- 
eral, the binary system cannot be in equilibrium in GR 
since GW carry angular momentum and energy away to 
infinity. However, the quasi-equilibrium approximation 
is well fulfilled since the time scale of evolution due to 
GW emission is longer than that of the orbital period 
even just before coalescence 0. In this situation, we can 
expect that (A) there exists an approximate Killing vec- 
tor £ expressing this symmetry of the space-time and the 
fiuid. Moreover, it has been pointed out by Kochanek 
and by Bildsten & Cutler [|| that (B) the flow fields of 
binary stars just before coalescence become irrotational, 
because (i) the viscosity of the neutron star matter is 
too weak to synchronize the spins and the orbital an- 
gular velocity by the time of merging and (ii) the initial 
spin angular velocity of each component star is negligible 
compared to the final orbital angular velocity. 

To construct realistic BNS's in quasi-equilibrium, the 
above two properties (A) and (B) should be implemented 
in the formulation of the basic equations. Strictly speak- 
ing, it is impossible to assume the property (A) in exact 
GR, because there arises GW emission from the source 
to the asymptotically fiat region. This implies that there 
is no unique way to formulate properly an approxima- 
tion which fulfills the property (A) for the Einstein equa- 
tions. However, several authors have developed formu- 
lations for BNS's in quasi-equilibrium. The formulation 
for the first and second post-Newtonian approximations 
(1PN/2PN, respectively) have been developed and co- 
rotating quasi-equilibrium configurations of BNS's in the 
IPN approximation have been calculated [Q. Also, the 
Wilson-Mathews formulation has been implemented by 
several authors to express quasi-equilibrium BNS's 

as well as space-time around them. In this paper we use 
the Wilson-Mathews formulation to express the gravi- 
tational field. For the fiuid, the property (B) is imple- 
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mented by introducing the velocity potential and assum- 
ing the barotropic flow. The property (A) is imposed 
by decomposing the variables into the direction of the 
Killing vector £ and the spatial direction pl|-p^. 

Consequently the basic equations we use in this paper 
are identical with those used by Bonazzola, Gourgoulhon 
and Marck jl^ and by Marronetti, Mathews and Wilson 
Ip^ . We will summarize the formulation and the basic 
equations below. Geometrical units with G = c = 1 are 
used throughout this paper. Latin indices run from 1 to 
3, whereas Greek indices run from to 3. 



B. Formulation for the gravitational field 

We adopt the ADM decomposition for the convenience 
of numerical computations. The metric is written as, 

= -^a^dt^ + 7y {dx' - uj'dt) {dx^ - uj^dt) , (1) 

where t/^^, a, tu' and 7^ are the 4-metric, the lapse func- 
tion, the shift vector and the 3-metric on a 3D spatial 
hypersurface S*, respectively. The unit normal to St is 
written as 

) and = (-a, 0,0,0), (2) 
a a J 



and is defined as 



(3) 



By using the (3+1) formalism, the Einstein equations 
are decomposed into constraint equations and evolution 
equations. The Hamiltonian constraint and the momen- 
tum constraints are, respectively. 



(4) 



DjK^^ - D'K = 8TTj\ (5) 
where the source terms pn and j'' are defined by 

PH = n^n'-'T^i,, (6) 



r 



(7) 



Here i?, Kij, K, Di and T^i/ are the scalar curvature, 
the extrinsic curvature, the trace part of Kij, the covari- 
ant derivative with respect to 7^, and the stress-energy 
tensor, respectively. 

According to the (3-fl) formalism, the 3-metric 7^ 
satisfies the dynamical equations 



d_ 
di 



7y = ~2aKij - DiUjj 



DjOJi. 



(8) 



These equations are decomposed into their trace part and 
trace free parts jl^ : 



-aK - DiUj' 



(9) 



-D.ujj - DjUj, + -^tjDkUJk- (10) 

where 7 = det7y and K = K^^. 

The dynamical equations for Kij are also derived and 
decomposed into their trace part and trace free parts. 
We only write the trace part : 



dK 

'dt 



aR - D'D.a + aK^ - uj'D.K 



'Aira (ipH — S) , 



where S is defined as 



s = Y'T,, 



(11) 



(12) 



The post- Newtonian approach up to the 2PN order ad- 
mits exact equilibrium states of BNS's |Q. However, the 
2PN equations for gravity require a significant amount of 
computation. Also, the problem of GR formulation for 
quasi-equilibrium states has not been settled. To make 
the problem tractable, we assume that the spatial part of 
the metric %j is conformally flat for all the time during 
the inspiraling stage as follows: 



-1/3 



lij — lij 



(13) 



(14) 



where /y is the flat space metric and ^ is a conformal 
factor §,0. This is a rather strong restriction for the 
space-time since it is well known that the 3-metric is not 
conformally fiat even for stationary axisymmetric space- 
time. However, under this assumption, we can compute 
exactly BNS's up to the IPN order and spherically sym- 
metric configurations in full GR. It is also expected that 
this choice minimizes the gravitational wave content of 
the (physical) space-time by removing the dynamical or 
"wave" degrees of freedom from the field. Moreover, 
since this choice can always be adopted to find initial 
data without any approximation, our solutions give ini- 
tial conditions for 3D simulations of BNS coalescence in 
fuU GR 0. 

For the time slicing condition we choose maximal slic- 



ing: 



= 0, 



(15) 



and we assume that this condition holds for all the time 
of the evolution 
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dK 
'dt 



0, 



(16) 



although the other evolution equations for Kij are omit- 
ted artificially. 

From the above assumptions Eqs. (|lj) - (|l6|), the basic 
equations for the gravitational field, become as follows. 
The scalar curvature reduces to 



R = 



(17) 



where = VV^ — P^ViV j is the Laplacian of flat 
3-space. Eqs. (||), (||), ( p^ ) and ( pT[ ) are rewritten as 



(18) 
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(19) 



if = Q (V'c^J" + V^c^') - ^f'Vk^'' ] , (20) 



(a*) = a* ( -'^-''K.jK'^ + 27r^'*(p + 25) ) . (21) 



Here we define K''^ as 



(22) 



and its indices are lowered by the flat metric fij , as Kij — 
fikfjiK^^ . Vi is the derivative with respect to the flat 
metric fij. 

To make Eq. (|l|) more tractable, the shift vector cj' is 
decomposed into the sum of a vector and the gradient of 
a scalar as 



1. 
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(23) 



Eq. (|l^) is divided into the following equations for these 
variables and B as, 

V^G* = -2Vj (a^--^) - 167ra*^i\ (24) 



(25) 



The basic variables of the gravitational potentials for 
the actual computation are then a, and B, and 
the corresponding equations are Eqs. (Oj}, (^), ( |2^ ) and 
(p5|), together with relations (20) and (|23|). It should be 
noted that all of these equations are elliptic type PDEs. 
We can solve these equations iteratively with appropriate 
boundary conditions by using a certain computational 
method, which is described in a later section. 



C. Formulation for the fluid part 

General relativistic irrotational flows have been con- 
sidered as the extension of the Newtonian case and ap- 
plied for accretion onto a black hole . The formu- 
lation of irrotational flow for inspiraling BNS's in quasi- 
equilibrium has been derived independently by Shibata 
Jin and by Teukolsky [Q. We will follow their formula- 
tion. 

The energy momentum tensor for a perfect fluid can 
be written as 

T^u = p{l + £ + u^Uu + Pg,iu, (26) 

where p, e, P and are the rest mass density, the spe- 
cific internal energy, the pressure and the 4-velocity, re- 
spectively. We assume a polytropic EOS for the fluid 



p = (r- i)p£ = Kpi+i/'^ 



(27) 



where n is the polytropic index, T — 1 -I- 1/n, and k is 
a constant. We define the relativistic specific enthalpy h 
as 



h=l + e+- = l + Khi+l) pi/". 
P 



(28) 



The conservation equation for the energy momentum 
tensor 

(4)v^T,^ = 0, (29) 

can be rewritten as Eulcr's equation 

("i^V^ (hu^) + (-^'v^/i = 0, (30) 

where ''^■'V^ is the covariant derivative with respect to 
5^^. The conservation of rest mass density is written as 

(4)V^ (puf^) = 0. (31) 
The covariant definition of irrotationality in 4D is 



Ufi 



= i ( W V. {hu^) - (4)v^^ ^ 0, (32) 

where P'^ = g^^ + u^Ui, is the projection tensor and we 
have made use of Eq. (^0|). Thus the quantity hU)_i can 
be expressed in terms of the gradient of a potential $ as 



(33) 



In the following, we refer to <& as the velocity poten- 
tial. Then the equation of rest mass conservation (|3^) is 
rewritten as 

(4)VM (4)y^$ ^ „ (4)yM (^Ij, £j (4) V^$. (34) 
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From the normalization condition for the 4-velocity 



-1, we get the fohowing equation for h 



h= -(4)v^$(4)yM$ 



1/2 



(35) 



This is the generahzed BernouUi equation in GR since 
it is also derived by integrating Euler's equation ( ^o|) as 
follows, 



(36) 



Now we decompose Eqs. (|3J) and ( |35| ) into the ADM 
(3+1) form and impose the quasi-equilibrium condition. 
The quasi-equilibrium condition is implemented by as- 
suming the existence of a timclike Killing vector £^ such 
that 







and 



(37) 



(38) 



where / denotes a certain fluid variable such as p, h or 
u^. Ci denotes the Lie derivative along Since we 
have assumed a polytropic EOS and irrotational flow, 
the number of independent fluid variables is two, i.e. the 
velocity potential <i> and, for example, the density p. Note 
however that not $ but its gradient is the fluid variable. 
Therefore, the following relation holds 



that is, 



(39) 



(40) 



where C is a positive 'global' constant on the fluid. 

Performing the (3-1-1) decomposition of Eqs. ( ^4| ) and 
( p5| ) , and imposing the quasi-equilibrium conditions ( |3^ ) 

nations for irro- 



and <M\), we have the following basic ec 



tational binary systems, respectively |1 1 ,|1 



a 
A 



D, In 



\2 



ap 



Here, B"^ is a rotating shift vector defined as 



(41) 



(42) 



(43) 



where is the generator of rotations about the rotation 
axis of the binary. A is defined as 



A = C + BW,^. 



(44) 



Using the quasi-equilibrium conditions Eqs. (|38|) and 
(^), we obtain the following expressions for the normal 
derivatives of a scalar quantity / and the velocity poten- 
tial $, 



1 



a 



(45) 



and 



= -- {c + BW^'^>) = 



A 

a 



(46) 



The boundary condition for the velocity potential 
equation is derived as follows from the fact that the spa- 
tial fluid velocity on the stellar surface flows along the 
stellar surface on each 3D spatial hypersurface. We de- 



compose u*^ as 



(#' + V''). 



(47) 



and assume that V^'' is a spatial vector, n^V^'' = 0. By 
using V^, the boundary condition is expressed as 



VW,p\ , = 0. 

I surf 

From the definition for y , we obtain 
1 



(48) 



\ ah 

Therefore is rewritten as follows 
_ 1 

„2 



(49) 



a 

T 



(50) 



where Eqs. ([4 3D, (|46|) and ( |49D are used. From this ex- 
pression, the boundary condition (Bq) is written as 



D'^ - ) D,p 



0. 



(51) 



This is identical to the one derived from imposing the 
finiteness of the r.h.s. of Eq. ( ^H ) on the compact support 
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of a star which assures that the elliptic PDE for $ is well- 
defined on it. 

Finally, we impose the spatial conformal flatness 
Eq. and the maximal slicing condition Eq. (|l5|). By 
using the following identities for the derivative of a scalar 
function, 



(52) 



DiA = V,A and D'AD,B = f^'^VMVjB, (53) 

the basic equations and boundary conditions for the fluid 
part, namely, Eqs. (|l|), (||), and (|l]) are rewritten as 
follows, 

2 A 
V'V,<I> = V**Vi$ + '^^B'\'^^ 



and 



where 



surf 



A = C + B^Vj^, 



(56) 



(57) 



We also need another boundary condition to determine 
the stellar surface. This comes from the definition of the 
surface itself. 



P = 0. 



(58) 



Eqs. (|5J) and ( pq ) form a system of elliptic type PDEs 
with Neumann boundary conditions. 

Next we derive the matter source terms which appear 
in the Einstein equations, i.e. pH^ j" and S. From the 
definitions, 



Ph = n^^n^T^^ 



p\^ 



P, 



S = Y'^T^. ^ j^-ph + 3P, 



(59) 
(60) 



/ = -f^n^T^,, ^ ^D'<S> = ■^^-^V*^. (61) 
ha ha 



III. THE COMPUTATIONAL METHOD FOR 
IRROTATIONAL BINARY SYSTEMS IN GR 

In this section, we describe the solution method for the 
basic equations derived in the previous section. This is 
essentially the extension of the method developed pre- 
viously by the authors and co-workers and successfully 
applied to solving for the structure of a rapidly rotating 
single star, or a binary system in Newtonian gravity or 
in GR 



eg = 7V/2 




sym 


Q; 


sym 


Lr 


sym 


G" 


ant i- sym 




sym 


sin if> 


sym 
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sym 



sym 
sym 



sym 1 + O 



sym l + O 



sym 



sym 



O 



sym anti-sym anti-sym O I — 



(54) 




df = n/2 = 




r-f = Rs 




P 


sym sym 


sym 


Eq. (g) 


(55) 


$ 


sym anti-sym 


anti-sym 


Eq. (g) 



TABLE L Symmetries of variables and their boundary con- 
ditions are listed : 'sym' and 'anti-sym' denote plane sym- 
metry and anti-symmetry, respectively. Boundary conditions 
are at the outer boundaries, i.e. Vg ^ oo for the gravitational 
potentials and rf = Rs for the fluid, where Rs is the stellar 
surface. The equatorial plane is taken to be the 9 = 71/2 plane 
and the stars are aligned along the {Og, ipg) — {tt/2, 0)-axis. 



A. Coordinates and symmetries 

For the present computation, we prepare two spheri- 
cal coordinate systems {r,9,ip), one for the gravitational 
field whose origin is at the intersection of the rotational 
axis and the equatorial plane, and the other for the fluid 
variables whose origin is at the coordinate center of the 
star. We will call them the gravitational coordinates and 
the fluid coordinates, respectively. When necessary, we 
will distinguish these two coordinates by subscripts g and 
/ as {rg,dg,ipg) 0,11(1 {v f , 9 f , (fi f) , respectively. 

We assume that the symmetry of the physical quan- 
tities is as shown in Table |. Note that the equatorial 
plane of the star is the 9g = 9f = tt/2 plane and that 
the axis of orbital motion is located at 9g = 0. Planes 
of (anti-)synimetry are the equatorial plane and a plane 
with ipg = for the gravitational coordinates, iff = 
and TT for the fluid coordinates. In this paper, we only 
consider an equal mass BNS system and hence a plane of 
(anti-) symmetry with respect io Lpg = it / 2 exists for the 
gravitational coordinates as well. 

Because of these symmetries, we only need to treat one 
octant of the spherical coordinate system for the gravita- 
tion, namely, {rg,9g, (pg) € [0, oo) x [0, 7r/2] x [0, 7r/2]. For 
the fluid coordinate system, we need to treat a quarter 
of the whole star, that is, (r/, 0/, tp/) € [0, i?] x [0, 7r/2] x 
[0, tt], where R is the largest coordinate radius of the star 
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in the fluid spherical coordinates. We set the coordinate 
center of the star to be at (r^, Og,ipg) = (d, n/2, 0) of the 
gravitational coordinate system, where 



d — (i?out + ^in)/2, 



(62) 



is defined as the coordinate center of the star, and 
{rg,6g,tpg) = (i^Qut , tt/ 2 , ) aud (i?in, 7r/2, 0) are the 
outer and inner edges of the star, respectively. We show 
a schematic figure of the coordinate systems in Figure |^. 

Since we have assumed conformal flatness in space, we 
may introduce these spherical coordinate systems in the 
same manner as we do for Newtonian models. Especially, 
since the derivatives appearing in the basic equations are 
those for the flat 3-metric fij, it is convenient to take 
the non-coordinate basis whose metric becomes the unit 
matrix fij — 6ij — diag( 1,1,1) as is usually done for 
vector analysis in flat space, where Sij is the Kronecker 
delta. We write the orthonormal non-coordinate basis 
for the spherical coordinates as 



{er,ee,e^}, where 



(63) 



When necessary, we will distinguish the bases for the 
gravitational coordinate system and the fluid coordinate 
system by the superscripts g and /, respectively. Using 
these bases, the gradient of a scalar function cf) is written 
as 



or 



Idcj) 



1 d(j) 
V sin0 dip 



(64) 



Hereafter, we often for convenience make use of the index 
notation and the notation of vector analysis in the same 
equation, as above. 



B. Poisson solver 

Seven elliptic type PDEs appear in our formulation. 
They are Eqs. (|l|), iM), (||) and (|2|) for the gravita- 
tional field, and Eq. ( p4| ) for the fluid variable. To solve 
these equations, the flat part of 3D Laplacian is first sep- 
arated out. Then, the equations are transformed into 
integral forms by using the Green's function of the 3D 
Laplacian. 



(a)- 



(65) 



Here the suffix (a) is an index to specify each equation, 
and A is the 3D flat Laplacian in spherical coordinates. 



A(/)(q) 



1 d 



1 



d 



• sin 6 do 



do 



1 



r2 sin dip^ 



(66) 



Coordinates for the 
gravitational field 



Coordinates 
for the star 





R g.mid 



Equatorial plane 



(9 = Til2) 



FIG. 1. Schematic figure of the coordinate systems used in 
the actual computation, (a) A sketch of the coordinate sys- 
tems for the star and the gravitational field, (b) The larger 
spherical coordinate system corresponds to that for gravita- 
tion (white and dark gray region) . The smaller spherical coor- 
dinate system corresponds to that for the fluid (gray region). 
Ro is defined as Ro = {Rout — Rin)/2. 
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By using a Green's function which satisfies the following 
equation 



Eq. (pq) is transformed as 



0(a) 



An J |r-r'| 



(67) 



(68) 



where x is either a function or a constant depending on 
the boundary condition. We use the Legendre expansion 
to compute the integral in Eq. (pq), 



1 



£/„(r,r')X^6„ 



(n — m)! 
(n + m)! 

xP„"(cos6') F„'"(cos6i') cosm(^ - ip'), 
where fn('r,r') is defined as, 



(69) 



Ury) 



£„,, is defined as, 



1 r' 
r \ r 



-f-V 

r' \r'J 



for r' < r , 
for r < r' , 



(70) 



for m = , 
for m = 1,2,. 



(71) 



and P„™(cos(?) is the associated Legendre function. Here 
dV = r'2 sin 6' dr'de' dip'. 

Eq. ( |6^ ) is a formal solution in a sense that each source 
term Sj^a) ^^^o contains the variables {(j)^a}}- Therefore, 
an iterative method is used to find a solution. All of 
the equations are consistently solved by a self-consistent- 
field (SCF) iteration scheme poj-p^, whose procedure is 
described in detail in Paper I of H 



C. Solution method for the gravitational field 

To apply the above procedure, we need to separate out 
the 3D flat Laplacian A for the spherical coordinates as 
in Eq. (|66|). We may regard operating on scalar val- 
ued functions, as in Eqs. (^8|), ( pl| ) and (p5|), in the same 
way as A in Eq. (|6^). On the other hand, the Laplacian 
operates on a vector valued function on the l.h.s. of 
Eq. ( P4|) . Although it is possible to use the Green's func- 
tion of the Laplacian for a vector valued function, in- 
stead, we explicitly write down its components and de- 
rive the form of Eq. (|6^) in each expression, so that we 
can follow the same procedure as for the scalars. 

In the gravitational coordinates, this procedure goes 
as follows. 



V^G^' = V'V.G^' = A (Ces + G^e^ + G^e^) , (72) 



[V^G'Y = e, • A (G''e9 + G'^e^ + G'^e 



AG'' - ^ G 



1 dQv dG'^ 



G« 



tan 6* sin6' dp 89 



(73) 



[V^G^] ' = ee • A (G'-e^ + G^e^ + G'^eS ) 



AG*^ 



G^ 



cos 61 dG'P dG' 



2 sin^ e sin"' 



dip 



(74) 



[V^G^] = • A (G'-eS -I- G^e^ + G^^e^ 



AGf - — 



G'^ 



1 dG'' cose dG^ 



\ 2 sin^ 9 sin 6 dp s\v? 6 dp) 



, (75) 



where Eqs. (^ - (|7|) are the l.h.s. of Eq. (H). For 
the p component of Eq. ( |24|), we multiply by sin 95 and 
rearrange the terms of Eq. ([75|) as follows: 

sin p = sin ^ • A (G'^eS + G^e^ + G^^e^ ) 



A(sin^G'^) 



sin^ ( 



cosf/smiyS- 



smf^sm^s- 



dp 



dG^ 
dp 



cos p- 



dGf 



dp 



(76) 



In these expressions we suppress the subscript g for the 
gravitational coordinates. 

Terms appearing in expressions ([73|), ( [74[ ) and (jT^) in 
addition to the flat 3D Laplacian A are transferred to the 
r.h.s. and are included as a part of the source term of each 
component. By the above rearrangement, all of the basic 
equations for the gravitational part have been written 
in the form of Eq. (|65|). Therefo re, w e can follow the 
procedure described in subsection III E| . As we show in 
Table ||, the boundary conditions arc imposed by setting 
X = 1 for Eqs. (|l|) and and x = for Eqs. (|2|) 
and (|2^). Because of the symmetries for the variables 
shown in Table |, the source terms have corresponding 
symmetries as well. As a result, appropriate behaviors 
of integral terms in Eq. (|6^) at r ^ 00 are automatically 
satisfied. 



D. Solving method for the fluid part 

The basic equations and boundary conditions for the 
fiuid part are Eqs. (|5^) - (|5|). Since these equations 
are essentially the same as those of the Newtonian case 
solved in Paper I, we can apply all of the computational 
techniques developed in that paper. As we mentioned 
in subsection [II A , we prepare the spherical coordinate 



system to compute the structure of the star whose origin 
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is at the coordinate center of the star. The elliptic type 
PDE Eq. is solved in this coordinat e syst em by us- 
ing the same procedure as in subsection IIIB. Eq. (jsj) 
has the same form as Eq. ( |65|) and is integrated to give 
Eq. (IH) . To impose the boundary condition Eq. (js^) , we 
regard the term x in Eq. ( |68|) as a regular homogeneous 
solution of the Laplace equation inside the star, i.e., 



For the numerical computation of the fluid part, it is 
convenient to introduce the surface fitted spherical coor- 
dinate system as used in the Newtonian case The 
surface of one NS in the binary system can be expressed 
by a function Rs {9f,iff) even when the deformation of 
the shape is relatively large. By using this function, the 
new coordinates {rj, 0*pip*j) are defined as follows, 



(77) 



In the spherical coordinates, such a homogeneous solu- 
tion to Eq. ( IttI ) can be expressed by using the associated 
Legendre functions as follows: 

oo I 

;=i m=l 



(2; + i)(;-m )!y/' 

4tt{1 + m)l 



Pi {cos 0f) sin rrnpf , (78) 



where ajm's are certain constants. Here we have taken 
into account the regularity of the solutions at r^^ =0 and 
the symmetries listed in Table |[ The coefficients a;,„ are 
computed so that the velocity potential <f> satisfies the 
boundary condition ([S^). 

Eq. ( p5| ) is used to compute the density distribution 
within the star. For numerical computations, we intro- 
duce a function q which is similar to the Emdcn function 
of Newtonian polytropes defined as follows: 



P 

q = —■ 
p 

In terms of this function, we can express 



P = K-"q'^+\ 



1 + (n+ l)g. 



(79) 

(80) 
(81) 
(82) 



Substituting the last three expressions into Eq. (pq), we 
find 



'/ 



Rs{0f,ipf) 



, 0}=ef, and ip}^ipf. (87) 



The stellar interior with the assumed symmetry is 
mapped into the region {rj,9j,(p*^) g [0,1] x [0, 7r/2] x 
[0,7r]. The surface fitted coordinate is advantageous for 
computing numerical derivatives and imposing boundary 
conditions at the stellar surface. Implementation of this 
coordinate system is described in detail in Paper I. 



E. Normalization of quantities and choice of 
parameters 

Non-dimensionalization of variables and proper choices 
of parameters are important for the iteration scheme 
to obtain converged equilibrium configurations stably. 
For convenience in making numerical computations, we 
rescale the coordinate length so that the coordinate ra- 
dius of the star at the intersection of the surface and the 
coordinate line with {9g,ipg) = (7r/2,0) is unity, namely, 
we introduce a parameter Rq so that Rs{tt/2,0)/ Rq = 
1 = i?s(7r/2, 7r)/i?o in the fluid coordinates. 

By using this i?0: the basic equations for the gravita- 
tional potentials, Eqs. (|l|), (|l|), (||) and (||), whose 
forms are essentially the same as Eq. ( |6^ ) or Eq. (|6^), 
become as follows, 



</'(a) = X - 



1 

47r 



I ry 2 Q m 
^0 -^(a) 



dV' 



(89) 



1 



1 



A2 



1/2 



(83) 



Accordingly, the matter source terms Eqs. (p9[)-(|6 
are rewritten as follows. 





(1 + (n -I- 1) g) a2 



(84) 



where A 



Here S,"l stands for the non- 

(a) 



dimensionalized part of the matter source terms corre- 
sponding to Eqs. (|5|)-(|l|) or (||)-(|6|) and S^^^ denotes 
remaining terms. We explicitly show the dependence of 
these quantities on the parameters k and i?o- In the equa- 
tions for the fluid part, k and i?o are canceled out after 
rescaling and non-dimensionalization. By this rescaling, 
the coordinates and variables are transformed as follows, 



PH + 2S ^ K-"q" 



3A2 



{l + {n + l)q) a2 



A 



J = K q 



(2+(2n-3)<7) , (85) 



(l + (n + l)g), 



(86) 



= K.jRo, B = ^, ^ = ^: ^Ro- (90) 
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We may consider that k is used for non - dimension- 
alization of the variables, and Rq for rescahng of the nu- 
merical computations. However, since these two param- 
eters appear in the equations only through the following 
combination: 



Rq 



Rn 



(91) 



we can regard them as a single parameter when we com- 
pute each solution. In this paper, we only consider 
K — const sequences, which could be appropriate for in- 
vestigation of the final inspiraling stage of BNS's. (When 
constructing solutions, we could mimic physical effects 
such as heating by changing the value of k.) 

We then have three parameters Rq, and C in the 
basic equations. We need to impose three more condi- 
tions to specify them. For this purpose, we choose three 
locations where q has definite values. Namely, we set 
q — At the intersections of the surface and the co- 
ordinate line with (9g,ipg) = (7r/2,0), i.e. at the two 
points (ry,0y,(/7y) = (l,7r/2,0) and (l,7r/2,7r) and also 
set g = at a grid point where q takes the largest 
value inside the star ||23|] . Since we have introduced the 
surface Rs{9*f,f*f) in the surface fitted coordinates, the 
former two conditions are explicitly imposed by setting 
i?s(7r/2,0) = 1 = i?s(7r/2,^). 

The above three conditions are applied to Eq. ( ^ ) at 
three points to get equations for Rq, D, and C. It is known 
||l|,[0| that a and * are scaled as 



By using this rescaling, Eq. ( pq ) is written as. 



j2«0{l + („+l)g}2+(S$) " 



(93) 



We use three equations which are derived by imposing the 
above three conditions at three points in the star, namely, 
at the inner and outer edges and at the coordinate center 
of the star. These three non-linear algebraic equations 
with respect to Rq, f2 and C are solved by using the 
Newton-Raphson method and their values are updated 
through the iteration procedure. 

The above choice is known to make the iteration con- 
verge stably. This choice of parameters and the computa- 
tion scheme for them are essentially the same as those for 
rotating single stars [|l| or those for co-rotating BNS's 
M in GR. 



F. Discretization and numerical computation 

Since we have rescaled the length by Rq, the compu- 
tational domain is measured by this unit. For the gravi- 
tational coordinate system, the whole computational do- 



main is taken as {rg,9g,ipg) G [0,R, 



■5, max J 



[0,^2] 



[0, 7r/2]. It is discretized equidistantly for the 9g and (pg 
directions. For the fg direction, we divide the region 
into two parts as [0, -Rg,mid] (the white region in Figure 
1(b)) and [i?g,mid, -Kg, max] (the dark gray region in Figure 
1(b)). We discretize the region [0, i?g^inid] equidistantly, 
and the region [i?g, mid, ^g, max] non-equidistantly. By de- 
noting grid points as iri,9j,ipk), in which we suppress 
the subscript g, the discretization becomes as follows. 



TV .'■ 

g,mid 



where 



< < i?3,mid, (0 < J < iv;,„id) 

Ti+i - Ti = 5(ri - fi„i), where 



(94) 



i?<,,mid < h < i?s,max, (iV;„,id < J < ^^max) , (95) 



tt/2 

'no' 
a 



where 



0<9,<-, {0<j<N'g] 



7r/2 



where 



< < -, (0 < J < N^) 



(96) 



(97) 



(^92) In Eq. (95), ^(> 1) is a certain constant. 

The surface fitted fluid coordinate system is discretized 
equidistantly in all directions as follows. 



' j-i-i 



1 

TV? 



where 



< f* < 1, (0 < i < iv;) 



7r/2 



where 



o<e* < -, (o<j<iv;), 



where 



< ^* < ^, (o < J < iv;) 



(98) 



(99) 



(100) 



The fluid coordinate system is always set within the re- 
gion [0, i?g,mid] of the gravitational coordinate system in 
order to maintain an accuracy. Standard finite differ- 
ences are applied for the basic equations with these grid 
points, which maintain at least second order accuracy. 

In the present computation, we typically took the val- 
ues listed in Table ||. A finer mesh is used for the region 
[0, i?g,mid] whose size in the Vg direction is about 2.5 times 
the diameter of one NS for the cases with types S, M, 
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g,mid 




5, mid 


' g.max 






Nf 


iV« 


NJ 


type 


5 


100 


40 


20 


20 


30 


8 


8 


16 


S 


5 


100 


60 


20 


30 


45 


12 


12 


24 


M 


5 


100 


80 


20 


40 


60 


16 


16 


32 


L 


4 


100 


80 


20 


40 


40 


20 


40 


28 


Ls 



TABLE II. Parameters for the computational domains and 
the numbers of grid points in each coordinate system are 
hsted. 



and L and 2 times for the case with type Ls. The region 
[Rg, mid, Rg, ma^] IS discietlzed Into 20 grid points in the 
•Tg direction. Types S, M, and L are mainly used to show 
the convergence of the scheme as functions of the mesh 
size and we will mention this in a later section. The mesh 
sizes of type S and type M in the region < < Rg.mid 
are, respectively, 1/2 and 3/4 of that of the finest type L. 
Type Ls is used for compute the configurations for rather 
small separations, although the numerical results of type 
L and type Ls are almost the same. 

Strictly speaking, Rg^max should be infinite in order to 
impose asymptotically flatness as a boundary condition. 
This could be implemented by using a certain appropri- 
ate coordinate transformation to compactify an infinite 
region, which has been used by several authors p^ , p5t . 
Although this is straightforward, we do not implement 
such a special treatment for simplification in the present 
numerical code. Instead, we truncate the computational 

It is known that such trunca- 



■g.mavi- 



domain at large R 
tion does not affect the local 
or the structure of each NS 



g,max 



)ropertics of the space-time 
!6[ . Although the choice of 
would affect the values of the phys- 



the value of R, 
ical quantities very little, it would be safer to cover the 
space as extensively as possible. We take a rather large 
value for Rg^max — 100 as shown in Table ||. Because of 
this choice, the error introduced into an integral qua ntity 
such as the ADM mass by this truncation (see Eq. (104) 
below) does not affect the digits of the numerical results 
presented in the later sections. 

We also truncate the order of the Legendre expansion 
in Eq. (^9|) and Eq. (|7^) at finite values instead of infin- 
ity. We denote the maximum value of the expansion for 
quantities in the gravitational coordinates in Eq . (|69| ) as 
n-max and that for the fluid coordinates in Eq. (|69|) and 
Eq. ( [78| ) (i.e. the expansion for the velocity potential) as 

^max- Typically, we use (^Trnaxi hnax 

) = (32,12). There- 
fore sources of truncation error are the finite difference 
and the truncation for this expansion. The results de- 
pend very little on l^^ax when we choose Zmax > 8. We 
will show how the results depend on rimax in a later sec- 
tion. 

As was mentioned before, the basic equations are 
solved iteratively. The iteration procedure is the same 
as that for Newtonian irrotational BNS models described 



in Paper I. The only difference in the numerical scheme 
is that an interpolation of physical quantities from one 
coordinate system to the other is required. We imple- 
mented the cubic spline interpolation method for all di- 
rections of the coordinates. As we will show in the 
later sections, one of the interesting problems relating 
to the irrotational binary sequence is to determine the 
behavior of the maximum density as the binary stars ap- 
proach. To investigate this problem, it is quite impor- 
tant to compute a smooth density distribution. In fact, 
an interpolation scheme of 1st order accuracy does not 
give a smooth density distribution, especially around the 
central region. Therefore, we have employed the cubic 
spline interpolation which keeps smoothness of interpo- 
lated data, although the cost for numerical computations 
is rather large. 

In our formulation, each quasi-equilibrium configura- 
tion is specified by three parameters, namely, the half 
of the separation between the coordinate centers of the 
two component stars d, the largest value of on the dis- 
cretized grid, and the polytropic index n. (Note that k 
has been scaled out.) The initial guess for the numer- 
ical iteration may be, for instance, an analytic solution 
of a Newtonian spherical polytrope with n = 1 at larger 
separation, say d = 2. We have started the iteration by 
setting Qc to be a small value and could then follow iter- 
ation cycles. We have regarded the obtained solution as 
being a converged one when relative differences of physi- 
cal variables of subsequent iterations on each coordinate 
grid point are less than a certain small value, typically 
we choose 10^'"'. Once we get a converged solution, we 
can change the values of parameters by a few tenths of 
a percent. Typically, a solution is converged after 100 
~ 300 iteration cycles. It takes about 1 minute for one 
iteration cycle using a DEC Alpha-station 1/166. 



IV. NUMERICAL RESULTS FOR 
IRROTATIONAL BINARY SYSTEMS IN GR 

A. Integrals 

In this section we show several integral quantities of the 
equilibrium configurations. The total rest-mass energy 
Mo, tot of the two component stars is defined as 

Mo,tot= / = / pu^" i-rif,^) dV, (101) 

Jm Jm 

where M. denotes the integration region which is the sup- 
port of the two component stars. The non-dimensional 
form of the total rest-mass can be written using Eq. ( |90| ) 
and (^) as 

Mo.tot = ^""/'Mo.tot = K-^'^Rl I K-^q^^^^dV 

Jm na 



^3 r g"A^'6^^ 
° Jm ha 



(102) 
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where dV = dV/ . The total mass-energy (ADM mass) 
is 



Ma 



DM 



._L / \/'^lfd^S^ V^^dK 



2tt 



2n 



= ^-^K,,K'^dV+ [ ^^pHdV (103) 

iDTl" Joe 

where Eq. ( p^ ) derived from the Hamiltonian constraint 
is used. This can be rewritten in the non-dimensional 
form as, 



Madm = «~"/'Madm = ^ / ^^'K,,.K''dV 



Ro 
16^ 



A2 



(104) 



The angular momentum is aligned with the dg — axis 
and can be defined as (see, e.g., ||27[|) 



tot — — f fijCK' d Sk o 

TT ./no O 



/ '^^"hCfdV^ f ^^"rgSinOgj'^dV, (105) 

JA4 JM 



where we used VkK^^ = '^^°DkK^'^ as well as the mo- 
mentum constraints (^, and £, — rg sin^^e^ from the def- 
inition. This is the total angular momentum contained 
in the space-time including both the orbital and spin an- 
gular momentum of the stars. Finally, we can substitute 
Eq. (|8^) for and write the angular momentum in the 
non-dimensional form 



Jtot = K"" Jtot = Rt / Sin eg3'^ dV 

J M 

= Rt ^'V„sin6l„^V'^$dF. 
Jm ha 



(106) 



In the following we will denote a half of the total rest- 
mass, the ADM mass and the angular momentum by 
Mo = Mo,tot/2, M = Madm/2 and J = Jtot/2, re_spec- 
tively. In the limit of large separation, Mq and M ap- 
proach the corresponding values for isolated stars. 

We also use the averaged separation da which coin- 
cides with the separation between mass centers of the 
two component stars in Newtonian limit, namely, 



da = 

Ro 



Mo 



r„ smt^o cos ipo — ; aV . 

ha 



M 



(107) 



B. Test of our new computational method 

In this paper we will show results for n = 1 polytropes. 
For this we need to specify two parameters in order to 



obtain one quasi-equilibrium configuration, i.e. the half 
of separation d and the largest value of g = gc- In or- 
der to check how our new computational code works, we 
have compared our results with (a) those of co-rotating 
GR binary systems tabulated in (lO| , and (b) those of the 
Newtonian irrotational binary systems in Paper I. With 
the former comparison (a), we can check the gravitational 
field part of our code for the strong gravity regime. We 
have also performed (c) a convergence test for the New- 
tonian case since it is desirable to show how the scheme 
actually converges to some analytic or semi-analytic re- 
sults. 

Concerning irrotational binary systems, Taniguchi p^ ] 
has analytically solved semi-detached irrotational binary 
systems with the IPN approximation for n — Q poly- 
tropes and Taniguchi and Nakamura recently developed 
a method to treat the compressible Newtonian case [ p9[ . 
Lombardi, Rasio & Shapiro have calculated the irrota- 
tional polytropic binary systems by partly including IPN 
correction terms |30[| . It will be important to perform 
the convergence test with these results as well as making 
a comparison with numerical results obtained by other 
authors who have used the identical formulation but dif- 
ferent numerical computational methods. 

It is convenient to characterize solutions or sequences 
by using the total rest mass Mo = const or, equivalently, 
the compactness {M/R)oo of an isolated spherical star. 
In actual computations, we choose a certain value for Mq 
and adjust the parameter qc during the iteration so that 
the value of Mo is kept constant. 



1. Comparison with other works 

Below, we show the results of comparison tests (a) 
and (b). In this section, type Ls in Table || is used for 
the numerical grids. For the expansions in Eq. (|69[) and 



Eq. (|78|), we include terms up to (nmax, ^max) = (20, 12). 



In Table III, we compare our results for co-rotating 
binary systems with those in ||l^. For these models we 
have replaced the equations for the fluid part by those 
for co-rotating fluids. Therefore, the code for the grav- 
itational part has been checked. We set the values of d 
and qc to be the same as those in [|lOf for each model. 



As shown in Table III, the calculated quantities VL and 



Ro which characterize the models for BNS's are in good 
agreement (relative differences are less than 1%). For the 
integral quantities Mqj M and J, the relative differences 
are at most 2%. 

Next, we compare our results for a quasi-equilibrium 
sequence of irrotational binary systems in weak gravity 
with those of Newtonian computations. It is reasonable 
to assume that the total rest mass is constant throughout 
the final evolution. Thus the quasi-equilibrium sequence 
with Mo = const mimics the realistic evolution sequence 
as a result of GW emission. In Figure 0, we compare our 
present results for the compactness [M/R)oo — 0.001 
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d 


n 

He 


Q 






M 


J 


0.20 


1.50 


0.0685 


0.079 


1.004 


0.1118 


0.10552 


0.02729 








0.079 


1.005 


0.1125 


10619 


02780 


0.00 


1.00 


0.0658 


0.101 


1.177 


0.1118 


0.10551 


0.02715 








0.101 


1.176 


0.1128 


0.10637 


0.02760 


0.20 


1.50 


0.2450 


0.168 


0.646 


0.1781 


0.16013 


0.04948 








0.168 


0.643 


0.1801 


0.16169 


0.05050 


0.00 


1.00 


0.2164 


0.202 


0.780 


0.1781 


0.16018 


0.05024 








0.203 


0.774 


0.1807 


0.16227 


0.05132 



TABLE III. Comparison between the results of Baumgarte 
et al. (1997) and those obtained with the present method for 
co-rotating GR binary systems. The first hne of each model 
corresponds to the results of Baumgarte et al. (1997) and 
the second line to the present results. The first two models 
correspond to the system with {M/R)oo =0.1 and the latter 
two models to {M/R)oo = 0.2 in Baumgarte et al. (1997). 
za = 0.0 corresponds to a contact binary system and za = 0.2 
to a semi-detached binary system where za = -Rin/^out . 



< 




6.0e-03 7.0e-03 „,^'0e-03 

Angular velocity 



with those in Paper I. In this figure, J is plotted against 
the separation dc or Cl. Relative differences of J at a 
fixed value of da or are ~ 0.5%, relative differences of 
da at a fixed value of J are ~ 1%, and relative differences 
of fj at a fixed value of J are ~ 1.7%. Therefore, our new 
code has reproduced the same results as those of other 
computations with reasonable accuracy. 



2. Convergence test 



3.2e-05 



c 

E 3.1e-05 

o 

E 



JO 
Ul 

< 



2.9e-05 




(M/R)_=0.001 
Newtonian 



3.8 4.0 _^,|.2 

Separation d^K " 



In order to see how well our numerical solutions ap- 
proximate true solutions, we have estimated the influence 
of the mesh size on the results. This can be done by in- 
vestigating how the values of physical quantities behave 
when the mesh size is varied. For this convergence test, 
we have used types S, M, and L in Table || for numer- 
ical computations. We chose (nmax, ^max) = (20, 12) for 
the expansions. In Fig. H, we plot the relative differences 
of the angular momentum J and the angular velocity 
51 from those of the semi-analytic results of the ellip- 
soidal approximation pl[ | for types S, M and L with val- 
ues of the separation dJUo being fixed. The plots clearly 
show that for global quantities, our results tend to con- 
verge as the mesh size is decreased. Note that even for 
d/Ro — 2.0, the difference between the polar radius and 
the equatorial radius is ~ 3% and this means that the 
results of the ellipsoidal approximation are reasonably 
accurate. 

On the other hand we need to be more careful about 
the convergence of local quantities such as the density 
distribution, the shape of the star, or the metric poten- 
tial distribution. Since we will discuss dependence of the 



FIG. 2. Comparison between the results of Paper I (New- 
tonian) and those of the present computation (with compact- 
ness {M/R)ao = 0.001). Quasi-equilibrium sequences for ir- 
rotational BNS systems with weak gravity (curve with dots) 
and with Newtonian gravity (curve with crosses) are plotted. 
The angular momentum as a function of the angular velocity 
(upper panel) and the angular momentum as a function of 
the separation (lower panel) are shown. 
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maximum density on the binary separation in a later sec- 
tion, we here discuss the property of the convergence of 
the maximum density as a representative of local quan- 
tities. 

In Fig. ^(a), we plot the maximum rest mass density 
Pmax against the half separation d/Ro, normalized by the 
Rq of each model, to show the convergence of solutions as 
the mesh type is changed. In the figure, the curves cor- 
respond to the mesh type S (the line with crosses), the 
mesh type M (the line with diamonds), and the mesh 
type L (the line with filled circles), respectively. We 
again set (nmax,^max) — (20,12). The maximum den- 
sity Pmaxid/Ro) does not change monotonically as the 
separation changes, and does not tend exactly to the 
asymptotic value /9max,oo, which is the central density of 
an isolated spherical star computed from the TOV equa- 
tion. However, the curves oscillate around Pmax.oo and 
the differences from it are at most 1% even for the small- 
est mesh type S. Note that the resolution of the star in 
the gravity coordinates becomes coarser for larger sep- 
arations. For instance, the mesh number of the gravity 
coordinates which cover the neutron star radius is only 6 
points in the ipg coordinate at d/Ro ~ 3.0 for the mesh 
type S. 

To show that pmax{d/Ro) converges to a certain num- 
ber for each separation, we plot the relative difference 
of Apmax/Pox against the mesh size, in Fig. ^(b). 

Here the relative difference is measured in terms of the 
value Pox extrapolated from the values obtained for the 
mesh types S, M and L in Table |l| for each d/ Rq. This 
figure shows that the local quantities such as the density 
also converge as the number of mesh points is increased. 

The relative difference of the maximum density is 
roughly proportional to (Ar)^ around d/ Ro ~ 1.5. This 
is satisfactory since we are interested in configurations 
with d/ Ro < 2.0, since the deformation of a star is small 
even at d/Ro ~ 2.0. For example the difference between 
the major and minor radii is only a few %. Therefore, an- 
alytic or semi-analytic approximations could be applied 
to configurations down to this separation |28|-30|. The 
results also indicate that the truncation error of the ex- 
pansion TT-max IS negligible for the case with d/Ro 1.5 
when we choose rimax ^ 20. 

To clarify the very small non-monotonic behavior of 

Pmax 

(d/Ro), we have computed sequences with AI/R — 
0.001 by changing the number of the Legendre functions 
in the expansion for the gravitational potential, rimax, 
in Eq. (p9|). In Figure ||(a), we plot Pmaxid/Ra) for 
?^max = 16, 20 and 24, using the mesh type M. For 
f^max = 16, the relative difference from the asymptotic 
value becomes ^ 1.5% at larger separations d/Ro ^ 2.5. 
On the other hand, the differences between the Umax = 20 
and n,„ax = 24 configurations are less than 0.5% ev- 
erywhere and Pma^id/ Ro) at larger separation almost 
reaches the asymptotic value Pmax.cx)- From this figure we 
may conclude that the number of the Legendre functions 
in the expansion, rimax, should be larger than 20. We also 
note that there arises a restriction for rimax and Zmax from 
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FIG. 3. Results of the convergence test for the weak grav- 
ity configurations (compactness {M/R)oo = 0.001). Rela- 
tive diflterences between the present results and those of the 
semi-analytic calculations with the ellipsoidal approximation 
by Lai, Rasio and Shapiro (1994) are plotted. The points on 
each line correspond to the results for the mesh types S, M 
and L in Table O. We plot the mesh size of the r-direction on 
the horizontal axis. The short dashed line is proportional to 
Ar, the long dashed line to (Ar)^, and the dash dotted line 
to {Arf. 



the number of grid points, since the number of numeri- 
cal grid points should be enough to resolve the periodic 
behavior as well as the orthogonality of the Legendre 
functions and trigonometric functions used in Eq. ( |69| ) 
and Eq. ([zs]). In Figure ||(b), we plot Pmax(d/i?o) for 
the n„iax = 20 and 32 sequences computed by using the 
mesh type L. The result with rtmax = 32 shows that the 
density is almost constant for 1.8 d/Ro ^ 2.2. We 
stress again that the present numerical computation is 
targeted to compute highly deformed configurations of 
component stars with d/Ro <, 2.0. Since all of the curves 
behave almost similarly for d/Ro <^ 2.0, our computa- 
tions are accurate enough to discuss the increase or de- 
crease of Pmax{d/Ro) of the solution sequence. Therefore 
a reasonable choice is rimax ^20. It should be noted 
that the integrated physical quantities depend little on 
the number rimax if "•max 20. 

Concerning the ?max, that is the number of terms in 
the expansion for the velocity potential <&, this does not 
affect local or global quantities when we choose ^max > 8. 
For instance, the choice of grid type L and Zmax = 12 is 
enough to compute $ accurately at any separation. 



C. Quasi-equilibrium sequences of irrotational 
binary systems in GR 

In Figures I and we show contours for GR irrota- 
tional BNS systems. We have used mesh types L and Ls 
for the computations in this subsection. As seen from 
these figures, our numerical method can handle highly 
deformed configurations of BNS systems and the grav- 
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FIG. 4. The convergence test for the maximum rest mass 

density Pmax for the case with compactness {M/R)ac ~ 0.001. 
(a) pina.x{d/ Ro) is plotted. Each line corresponds to a different 
mesh size, the mesh type S (the line with crosses), the mesh 
type M (the line with diamonds), and the mesh type L (the 
line with filled circles), respectively. The long dashed line 
shows pmax of an isolated spherical star pmax.oo- The short 
dashed lines show ±0.5% difference from this value, (b) The 
convergence of the relative error for pmax. The mesh size of 
the r-direction is used for the horizontal axis to represent the 
resolution. The short dashed line is proportional to Ar and 
the long dashed line to (Ar)^. 
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FIG. 5. The dependence of Pmax(d/i?o) on the order of 
Legendre expansion of the gravitational potential rimax for 
the case with {M/R)oo = 0.001. Each line corresponds to 
a different value of rimax. (a) The mesh type M is used for 
computations, (b) The mesh type L is used for computations. 
The horizontal long dashed line shows the pmax of an isolated 
spherical star pmax.oo. The horizontal short dashed lines show 
±0.5% difference from this value. 
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Equatorial plane 



FIG. 6. Contours of the conformal factor "^"^ (upper left 
panel) , the rest mass density distribution (upper right panel) 
and the shift vector with the stellar surfaces (lower right 
panel), for the model with (M/R)^ = 0.14 and d = 1.25 
in Xy-plane {9g = 9f — 7r/2-plane). Lengths in the figures 
are normalized with Rq. 



itational field around them. In particular, in the right 
panel of Fig. ^, we show the contours for the rest mass 
density p with {M / R)oo = 0.14. From this figure, we can 
see that a cusp- like structure appears at the inner edge of 
stars before the two stars make contact with each other 
p3| . The cusp is supposed to correspond to the inner 
Lagrange point (LI point), with mass overflow occurring 
from this point when the separation of two stars becomes 
smaller enough. 

In Table IV, we tabulate the values of physical quanti- 
ties for irrotational BNS systems with a cusp-like struc- 
ture for various degrees of compactness. We also plot 
the dimensionless total angular momentum ^tot/-A^ADM 
in Fig. H which is unity for an extreme Kerr black hole. 
For the co-rotating models [Q, Jtot/M'^jiyi becomes 
unity around {M/R)oc ~ 0.175. On the other hand, this 
value for irrotational BNS systems becomes unity around 
(M/i?)tx) ~ 0.12 as shown in Table IV. Since there exists 



a velocity field which rotates in the counter direction with 
respect to the orbital motion of each component star, the 
total angular momentum becomes smaller. Therefore, af- 
ter coalescence of the two stars due to GW emission, irro- 
tational BNS systems could form a single Kerr black hole 
by losing much less angular momentum from the system 
than co-rotating BNS systems. 

In Figures ^, we show quasi-equilibrium sequences 
for several values of {M/R)oo- The binding energy 
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FIG. 7. Structure of a component star of an irrotational 
BNS system. Contours of the distribution of rest mass density 
p (upper three panels) and of the velocity potential <I> (lower 
three panels) are drawn for the model with {M/R)oo = 0.17 
and d = 1.25. YZ-p\ane (upper left panel), ZX-p\ane (upper 
right) and XF-plane (lower right panel) correspond to the 
iff = and TT plane, the iff = 7r/2 plane, and the 6f = n/2 
plane, respectively. Lengths in the figures are normalized with 
Ro- 
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(M/R)^ 


Mo 


M 




Jtot /A/adm 


-Ro 


0.1 


0.112 


0.105 


0.0110 


1.08 


1.07 


0.12 


0.130 


0.121 


0.0148 


1.01 


0.979 


0.14 


0.146 


0.135 


0.0192 


0.963 


0.899 


0.17 


0.166 


0.150 


0.0271 


0.911 


0.783 



TABLE IV. Physical quantities of the irrotational BNS 
near to the configuration with a cusp-hke structure (LI point). 
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FIG. 9. The binding energy (upper panel) and the angu- 
lar momentum (lower panel) are plotted against the orbital 
angular velocity. Each curve corresponds to an Mo = const 
sequence. The value of the compactness {M/R)oo of an iso- 
lated star is attached to each curve. Cusp-like shape appears 
for the terminal model with the largest flAIo along each curve. 
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FIG. 8. Comparison of Jtot/Af^DM between co-rotational 
and irrotational BNS systems at the closest distances for var- 
ious values of the compactness parameter {M/R)oo- For the 
co-rotating models, values for the contact phases are plotted. 
For the irrotational models, those for the configurations with 
cusps are displayed. 



(M — Mod) /Mao and J are plotted against the normal- 
ized orbital angular velocity ^IMq. Since each curve is 
a sequence with constant rest mass Mq, an evolutionary 
track of a BNS system as a result of GW emission can 
be approximated by following each curve from smaller to 
larger values of ^Mq. The terminal points of the curves, 
i.e. the points with the largest ^IMq, roughly correspond 
to configurations with an LI point. In analogy with New- 
tonian models as well as with GR co-rotating binary sys- 
tems ||lo| , |3^ , turning points on these curves are expected 
to correspond to the points where dynamical instability 
sets in. Our results show that none of the systems with 
different values of the compactness {M/R)ao have turn- 
ing points. Therefore irrotational binary configurations 
with n = 1 are thought to be dynamically stable. In 
other words, at the final inspiraling phase, it is probable 
that the mass exchange starts earlier than the onset of 
the orbital motion instability. 

D. Individual collapse of the irrotational BNS 
system 

There are two main reasons why the close binary sys- 
tem of irrotational stars has attracted wide attention. 
One is that, as mentioned in Section & it is more proba- 



17 



ble that such a situation is reahzed just before coalescence 
due to GW emission, since viscosity is not strong enough 
to synchronize the spin of a NS and the orbital motion 
1^ at this stage. The other reason is that the numerical 
simulation of irrotational binary systems done by Wil- 
son, Mathews and Marronctti showed the result that 
the maximum density of each NS increases during the 
inspiraling due to the GR effect. They proposed a new 
scenario for the final inspiraling phase of irrotational BNS 
systems, namely, each component star of a BNS system 
individually collapses to form a binary black hole system, 
and after that two black holes would coalesce. 

This scenario has been debated by many authors 
(see |l^j3^ and references therein). Recent numeri- 
cal computations by Marronetti, Mathews and Wilson 
[ p^ again show the maximum density increasing by a 
few percent for a sequence of models with larger com- 
paction (M/R)oo = 0.19, although the most recent 
re-computations by Mathews and Wilson give a much 
smaller central density increase for polytropes ]36[ |. On 
the other hand, the computation by the Meudon numer- 



ical relativity group |15| does not show this tendency. In 
Figure 10, we show our results for the relative change of 
the maximum energy density (e,nax — eoo)/eoo of a star 
against the separation d. Here e = p{l + e), Cmax is its 
maximum value and e^o is its value in an isolated state 
which is computed by solving the TOV equation. 

We show the results for rimax = 20 in Fig. |lO|(a) and 
'^■max = 32 in Fig. p^(b). We notice that the results do 
not depend on the order of the Legendre expansion rimax 
for d/Ro < 2.0 with any values of the compactness. For 
d/Ro ^ 2.0, Cmax does not become constant for rimax = 
20 but tends to being almost constant when we choose 
«max = 32. To show the convergence around d/Ro <, 2.0 
clearly, differences between the result with rimax = 32 
and rimax = 20, (emax,32 - eniax,2o)/eoo, are plotted in 
Fig. 0(c). The relative difference (einax,32-emax,2o)/eoo 
is less than 0.35% at d/Ro < 2.0 for all (M/R)^ < 0.19. 

Differences of the relative change of the central den- 
sity (cniax — etx3)/eoo at larger separation d/Ro 2.0 get 
larger as the compactness becomes larger. For instance 
they are 0.5% and 1% for {M/R)^ = 0.17 and 0.19, re- 
spectively. However, in any case, differences are <^ 1% 
for any value of the compactness as shown in Fig. p^(b) 
and e„iax tends to a constant value when the component 
stars are detached. Therefore we may discuss the change 
of (cmax — eoo)/eoo along each solution sequence. Our re- 
sults show that the energy density changes in the range of 
0.5% to 2.5% and decreases as the separation decreases. 
Therefore, we do not observe the maximum density to in- 
crease and hence there is no tendency towards individual 
collapse in the present computation. Also the rate of de- 
crease of the maximum energy density around d/Ro ^5 2.0 
is larger for larger {A4/R)oo- 
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FIG. 10. The relative change of the maximum energy den- 
sity is plotted against the separation. The horizontal axis cor- 
responds to the coordinate separation in units of iio for each 
model, (a) The results computed setting the order of expan- 
sion to nmax = 20. Curves are plotted for (M/R)oo ~ 0.001 
to 0.19. (b) The results computed setting the order of expan- 
sion to Timax = 32. Curves are plotted for {M/R)oa ~ 0.001 
to 0.19. The maximum density does not increase as the sep- 
aration decreases, (c) The differences between the results 
with Timax = 32 and Wmax = 20. Curves are plotted for 
[M/R)oo ~ 0.001 to 0.19. These two results converge for 
d/Ro ^ 2.0. 
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V. DISCUSSION AND CONCLUSION 



results. 



A. Position of the present numerical method 



B. Future prospects 



As mentioned in section [I A, there are two different 



numerical schemes for solving irrotational binary star 
systems in GR [ p^p^ . The basic formulation and basic 
equations are almost the same but the numerical schemes 
are different. 

Bonazzola, Gourgoulhon and Marck ||l^ have used the 
multi-domain pseudo spectral method. They have in- 
troduced several coordinate systems one of which is a 
coordinate system fitted to the stellar surface. The basic 
equations which are essentially the same as those in this 
paper are solved by using the pseudo-spectral method on 
each of these domains. Although the implementation of 
this numerical method is rather complicated, it produces 
numerically 'exact' solutions with excellent accuracy. 

Marronetti, Mathews and Wilson |l^] have used a 
Cartesian coordinate system and a finite difference 
method. The equations which they solved are also the 
same as ours except for their approximate treatment on 
the boundary condition for the velocity potential equa- 
tion Although their computation is rather less ac- 
curate, Cartesian coordinates are commonly used in 3D 
numerical relativity to compute dynamical coalescence of 
compact binary systems, and hence it could be advanta- 
geous for providing initial data for future 3D simulations 
of BNS coalescence. 

Despite the existence of these two methods, we have 
presented another scheme in this paper. There are several 
reasons for this. Firstly, since quasi-equilibrium states of 
binary neutron star systems have an important signifi- 
cance in GW investigations, it would be better to get re- 
liable solutions by studying them from many directions 
even if this is the same problem. Secondly, results ob- 
tained from the two previous schemes seem to be some- 
what different. In particular, the behavior of the maxi- 
mum density during the evolution is different as discussed 
in subsection IV El| , although the recent re-computation 
by Mathews and Wilson [g6| gives consistent results with 
others for the polytropic configurations. 

In our solution method, we have used the integral form 
to solve Poisson type equations. A spherical coordinate 
system is introduced to implement this Poisson solver. 
From our experience with numerical computations, it 
seems desirable to use spherical coordinate systems to 
solve for the structures of component stars. On the other 
hand, to solve for the gravitational field, it is possible 
to implement Cartesian coordinates and suitable sophis- 
ticated Poisson solvers for the coordinates (for instance, 
the multi-grid method, ICCG, and so on). In our scheme, 
the coordinates for the star are patched on those for the 
gravitational field. This can be applied straightforwardly 
for a Cartesian grid. Such an implementation could also 
be useful for preparing initial conditions for numerical 
simulations as well as for checking the present numerical 



It is important to compute quasi-equilibrium configu- 
rations to connect the inspiraling phase, where the per- 
turbative expansion technique can be applied, to the 
merging phase, where fully dynamical computations are 
required. Such computations are also useful since the 
quasi-equilibrium solutions give accurate and appropri- 
ate initial conditions for simulations of BNS coalescence 

. Recently, a simulation of coalescing irrotational BNS 
systems in 3D full general relativity has been performed 
using the solutions of the present paper as initial config- 
urations ||l^ . The successful result of these simulations 
may be considered as showing that the present results 
are reliable for use as initial data for binary coalescence. 

In our formulation we have used the so called Wilson- 
Mathews formulation to treat the GR effects. However it 
is known that this approach is exact only up to IPN order 
for binary systems in terms of the post-Newtonian termi- 
nology, although 2PN or higher order corrections should 
be included to express realistic BNS systems. Asada and 
Shibata have formulated the Einstein equations so as 
to express BNS systems exactly up to 2PN order. Imple- 
menting their equations in numerical computations may 
be one of the possible ways to construct quasi-equilibrium 
BNS systems, although it is necessary to solve 29 ellip- 
tic PDEs simultaneously for the gravitational field. Re- 



3^ have developed a 
ibrium configurations 



cently, Usui, Uryu and Eriguchi 
new scheme to compute quasi-equif 
by using a quite different form of the metric from that of 
the Wilson-Mathews formulation. It would be interest- 
ing to compare results from that method with our present 
results. Also several authors have constructed a formula- 
tion for so called 'intermediate binary black hole' (IBBH) 
states This formulation could also be applicable to 
quasi-equilibrium states of BNS systems. 



C. Summary and conclusion 

We have constructed a numerical method to compute 
quasi-equilibrium states of irrotational binary systems 
composed of compressible gaseous stars in general rel- 
ativity. The new numerical method is based on a stan- 
dard finite difference method and therefore implemen- 
tation is rather simple and straightforward. We have 
used two coordinate systems for the computation, one 
for the gravitational field variables and the other for the 
fiuid variables. Physical quantities are interpolated from 
one to the other. We have shown that such a procedure 
works well by actual numerical computations of irrota- 
tional BNS systems in GR. 

The basic equation for the irrotational flow becomes 
an elliptic PDE with a Neumann type boundary condi- 
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tion at the stellar surface. For the treatment of general 
relativistic gravity, we have used the Wilson-Mathews 
formulation. In this formulation, the basic equations are 
again expressed by a system of elliptic PDEs. We have 
developed a method to solve these PDEs. We have cali- 
brated our solutions by using the results of independent 
numerical computations as well as of semi-analytic cal- 
culations. We have also tested the convergence of these 
numerical solutions. 

We have succeeded in computing quasi-equilibrium se- 
quences of irrotational BNS systems in GR. As in the 
Newtonian case, the BNS systems have been found to 
be dynamically stable during the inspiraling stage until 
Roche lobe overflow starts at the LI point. The quasi- 
equilibrium sequences also suggest that no substantial 
increase of the maximum energy density occurs during 
the inspiraling phase as a result of GW emission. Fur- 
ther investigations of numerical computations and im- 
plementation of various realistic EOS's and/or a realistic 
formulation of GR gravity are interesting and inevitable 
issues for theoretical predictions of GW signals supposed 
to be detected by the interferometric GW detectors in 
the next century. 
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